/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/****************************************************************************************************************
 * definition.h -- this module collects various constants, conversion factors and useful macro definitions. 
 ****************************************************************************************************************/

#ifndef DEFINITIONS_H
#define DEFINITIONS_H

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <signal.h>
#include <sys/time.h>

// Mathematical constants
#define	PI	3.14159265358979323846 // pi - fundamental constant

/* Physical constants */
#define	kB_IN_J_OVER_K		1.3806505e-23 /* [J / K] Boltzmann's constant */
#define	kB_IN_kJ_OVER_K		( kB_IN_J_OVER_K / 1000 ) /* [kJ / K] Boltzmann's constant */
#define	kB_IN_kJ_OVER_molK	( kB_IN_kJ_OVER_K * AVOGADRO ) /* [kJ / (mol*K)] Boltzmann's constant */
#define AVOGADRO		6.0221367e23 /* [mol^(-1)] Avogadro's number */
#define AMU		1.6605402e-27 /* [kg] Atomic Mass Unit */
#define EPSILON0_IN_F_OVER_m        8.854187817e-12 /* [F/m] = [C^2/(Nm^2)] = [C/(Vm)] vacuum permittivity */
#define DEBYE_IN_C_PER_m     ( 3.33564e-30 ) /* [Cm] - 1 Debye equals 3.33564*10^-30 coulomb*meter */
#define e_IN_C     ( 1.60217653e-19 ) /* [C] - "e" (charge of a proton) equals 1.602*10^-19 coulomb */

/****************************************************************************************************************
 * Brahms units (same convention used in the popular Gromacs program, see Gromacs user manual, ch. 2)
 * These constants must be used to convert any input parameter from physical units to internal Brahms units.
 ****************************************************************************************************************/
// fundamental units:
#define nm_IN_BRAHMS_UNITS      1. // length unit in Brahms is 1 nanometer (symbol: nm)
#define amu_IN_BRAHMS_UNITS     1. // mass unit in Brahms is 1 atomic mass unit (symbol: amu)   
#define ps_IN_BRAHMS_UNITS      1. // time unit in Brahms is 1 picosecond (symbol: ps)   
#define e_IN_BRAHMS_UNITS       1. // charge unit in Brahms is 1 electronic charge (symbol: e)   
#define K_IN_BRAHMS_UNITS       1. // temperature unit in Brahms is 1 electronic charge (symbol: K)   
// derived units:
#define fs_IN_BRAHMS_UNITS      ( 1e-3 * ps_IN_BRAHMS_UNITS ) 
#define kJ_PER_mol_IN_BRAHMS_UNITS       1. // energy unit in Brahms is 1 kJ/mol   
#define kcal_PER_mol_IN_BRAHMS_UNITS  ( cal_IN_J * kJ_PER_mol_IN_BRAHMS_UNITS )
#define J_IN_BRAHMS_UNITS         ( AVOGADRO / 1e3 )
#define C_IN_BRAHMS_UNITS   ( 1. / e_IN_C ) // 1 Coulomb in Brahms charge Units
#define Angstrom_IN_BRAHMS_UNITS  (1. / LENGTH_Angstrom )
#define m_IN_BRAHMS_UNITS   ( m_IN_nm * nm_IN_BRAHMS_UNITS ) //  THIS IS OK! 1 meter in Brahms length Units (1 m = 1e9 BU)
#define nm_IN_m            ( 1e-9 ) // 1 nm = 10^-9 m
#define m_IN_nm            ( 1. / nm_IN_m ) 
#define cm_IN_BRAHMS_UNITS   ( cm_IN_m * m_IN_BRAHMS_UNITS )
#define D_IN_BRAHMS_UNITS   ( DEBYE_IN_C_PER_m * C_IN_BRAHMS_UNITS * m_IN_BRAHMS_UNITS ) // 1 Debye in brahms units
#define DIPOLE_D                  ( 1. / D_IN_BRAHMS_UNITS ) 
#define D_IN_eA             (0.2082) // 1 D = 0.2082 electroUnit * Angstrom
#define s_IN_BRAHMS_UNITS   ( 1e12 ) // 1 second in Brahms time units 
#define bar_IN_BRAHMS_UNITS      ( 1e5 * N_IN_BRAHMS_UNITS / Sqr (m_IN_BRAHMS_UNITS) ) // 1 bar in Brahms pressure Units
#define atm_IN_BRAHMS_UNITS      ( atm_IN_bar * bar_IN_BRAHMS_UNITS ) // 1 bar in Brahms pressure Units
#define kg_IN_BRAHMS_UNITS           ( 1. / AMU_IN_kg ) // 1kg in brahms units, eg, in amu
#define g_IN_BRAHMS_UNITS           ( g_IN_kg * kg_IN_BRAHMS_UNITS ) 
#define N_IN_BRAHMS_UNITS            ( kg_IN_BRAHMS_UNITS * m_IN_BRAHMS_UNITS / Sqr( s_IN_BRAHMS_UNITS ) ) 
#define EPSILON0_IN_BRAHMS_UNITS  (EPSILON0_IN_F_OVER_m * (Sqr(C_IN_BRAHMS_UNITS) / (N_IN_BRAHMS_UNITS*Sqr(m_IN_BRAHMS_UNITS))))
#define SQRT_4_PI_EPSILON0_IN_BRAHMS_UNITS  ( sqrt( 4. * PI * EPSILON0_IN_BRAHMS_UNITS ) ) 

/****************************************************************************************************
 * The following constants are used to convert physical parameters/properties from internal Brahms 
 * units into physical units (for example when outputting properties of interest).
 ****************************************************************************************************/
#define LENGTH_Angstrom  (LENGTH_nm * nm_IN_Angstrom )
#define LENGTH_nm                 ( 1. / nm_IN_BRAHMS_UNITS )
#define LENGTH_IN_m       (LENGTH_Angstrom * 1e-10) // multiplier factor to trasform a brahms length in meters
#define BRAHMS_LENGTH_IN_cm ( 1./ cm_IN_BRAHMS_UNITS )
#define MASS_amu                  ( 1. / amu_IN_BRAHMS_UNITS )
#define BRAHMS_MASS_IN_g           ( 1. / g_IN_BRAHMS_UNITS ) 
#define TIME_ps                   ( 1. / ps_IN_BRAHMS_UNITS ) 
#define CHARGE_e                 ( 1. / e_IN_BRAHMS_UNITS )
#define CHARGE_C                 ( CHARGE_e * e_IN_C )
#define TEMPERATURE_K             ( 1. / K_IN_BRAHMS_UNITS ) 
#define BRAHMS_ENERGY_IN_kcal_PER_mol   ( J_IN_cal ) // assumes brahms energy unit is 1 kJ / mol  
#define kB_IN_BRAHMS_UNITS   ( kB_IN_J_OVER_K * J_IN_BRAHMS_UNITS ) // boltzmann's constant
#define PRESSURE_atm              ( 1. / atm_IN_BRAHMS_UNITS ) 
#define PRESSURE_bar              ( 1. / bar_IN_BRAHMS_UNITS ) 
 
// conversion factors
#define bar_IN_atm                0.986923267 // 1 bar = 0.986923267 atm
#define atm_IN_bar               ( 1. / bar_IN_atm )
#define kg_IN_g                   1000
#define g_IN_kg              ( 1. / kg_IN_g )
#define AMU_IN_kg        (1.660538782e-27) // 1 amu = 1.660538782e-27 kg
#define m_IN_cm                  100
#define cm_IN_m                  ( 1. / m_IN_cm ) 
#define nm_IN_Angstrom            10
#define cal_IN_J                 ( 4.184 )
#define J_IN_cal                 ( 1. / cal_IN_J )
#define m_in_cm                   100
#define VOLUME_nm3                (Cube( LENGTH_nm ))
#define VOLUME_m3                 (Cube( LENGTH_IN_m ))
#define AREA_nm2                  (Sqr( LENGTH_nm ))
#define AREA_m2                   (Sqr( LENGTH_IN_m ))
#define LENGTH_cm                 (LENGTH_IN_m * 1e2)  /* [cm] */
#define ENERGY_cal_mol            ( BRAHMS_ENERGY_IN_kcal_PER_mol * 1000 ) /* [cal/mol] */
#define ENERGY_J                  ( 1. / J_IN_BRAHMS_UNITS ) 
#define FORCE_kJ_mol_nm          ( BRAHMS_ENERGY_IN_kcal_PER_mol * cal_IN_J / LENGTH_nm )
#define FORCE_N                   ( ENERGY_J /  LENGTH_IN_m ) // Force [N]
#define FORCE_pN                  ( FORCE_N * 1e12 ) // Force [pN]
#define TIME_s                    ( 1e-12 ) // assumes ps is the fundamental unit of time in brahms
#define TIME_mus                  ( TIME_s * 1e6 ) // microsecond
#define TIME_ns                   ( TIME_s * 1e9 ) 
#define TIME_fs                   ( TIME_s * 1e15 ) 
#define DENSITY_rel               ( kg_IN_BRAHMS_UNITS * kg_IN_g / Cube( LENGTH_IN_m * m_in_cm ) ) /* relative density = density in g/cm^3*/
#define TENSION_dyn_cm            ( FORCE_N * 1e3 / LENGTH_IN_m ) // surface tension unit: dyn/cm = mN/m      - 1 dyn = 1e-5 N
//#define TENSION_dyn_cm            ( FORCE_N * 1e5 / LENGTH_CM ) // surface tension unit: dyn/cm = mN/m      - 1 dyn = 1e-5 N
                                             
/* miscellany */
#define MAX_SITES_PER_LIPID    (50) // Used in latDiff measurements, maybe it will need to be increased for very big lipids... 

//typedef float real;
typedef double real; // [Rapaport, p.23]
typedef struct { real x, y, z; } VecR; /* vector with real components [Rapaport, p.23] */
typedef struct { int x, y, z; } VecI;  /* vector with int components */
typedef struct { real x, y; } VecRTwo; /* two-dimensional vector with real components */
typedef struct { int i1, i2; } PairInt;  /* two-dimensional structure with two int components */

/* note: a matrix is treated as a linear array, e.g.:
   M = m[0] m[3] m[6]
   m[1] m[4] m[7]
   m[2] m[5] m[8] */

typedef struct {
  real u[ 9 ];
} RMat;

/****************************************************************************************************************
 * Macro definitions, mostly from Rapaport [p.481-488] 
 ****************************************************************************************************************/
#define Sqr( x )     ((x) * (x))
#define Cube( x )    ((x) * (x) * (x))
#define Sgn( x, y )  ( ( ( y ) >= 0 ) ? ( x ) : ( - ( x ) ) )
#define IsEven( x )  ( ( x ) & ~1 )
#define IsOdd( x )   ( ( x ) & 1 )
#define Nint( x )  ( ( ( x ) < 0. ) ? ( - ( int ) ( 0.5 - ( x ) ) ) : ( ( int ) ( 0.5 + ( x ) ) ) )
#define Min( x1, x2 )  ( ( ( x1 ) < ( x2 ) ) ? ( x1 ) : ( x2 ) )
#define Max(x1, x2)  (((x1) > (x2)) ? (x1) : (x2))
#define Min3(x1, x2, x3) (((x1) < (x2)) ? (((x1) < (x3)) ? (x1) : (x3)) : (((x2) < (x3)) ? (x2) : (x3)))
#define Max3(x1, x2, x3) (((x1) > (x2)) ? (((x1) > (x3)) ? (x1) : (x3)) : (((x2) > (x3)) ? (x2) : (x3)))
#define Clamp(x, lo, hi) (((x) >= (lo) && (x) <= (hi)) ? (x) : (((x) < (lo)) ? (lo) : (hi)))

/* basic vector operations: */
#define VSet(v, sx, sy, sz)  (v).x = sx, (v).y = sy, (v).z = sz
#define VSetTwo(v, sx, sy)  (v).x = sx,  (v).y = sy
#define VCopy(v1, v2)  (v1).x = (v2).x,	(v1).y = (v2).y, (v1).z = (v2).z
#define VCopyTwo(v1, v2)  (v1).x = (v2).x, (v1).y = (v2).y
#define VScale(v, s)  (v).x *= s, (v).y *= s, (v).z *= s
#define VScaleTwo(v, s)  (v).x *= s, (v).y *= s
#define VSCopy(v2, s1, v1)  (v2).x = (s1) * (v1).x, (v2).y = (s1) * (v1).y, (v2).z = (s1) * (v1).z
#define VSCopyTwo(v2, s1, v1)  (v2).x = (s1) * (v1).x, (v2).y = (s1) * (v1).y
#define VAdd(v1, v2, v3)  (v1).x = (v2).x + (v3).x, (v1).y = (v2).y + (v3).y, (v1).z = (v2).z + (v3).z
#define VAddTwo(v1, v2, v3)  (v1).x = (v2).x + (v3).x, (v1).y = (v2).y + (v3).y
#define VSub(v1, v2, v3)  (v1).x = (v2).x - (v3).x, (v1).y = (v2).y - (v3).y, (v1).z = (v2).z - (v3).z
#define VSubTwo(v1, v2, v3)  (v1).x = (v2).x - (v3).x, (v1).y = (v2).y - (v3).y
#define VMul(v1, v2, v3)  (v1).x = (v2).x * (v3).x, (v1).y = (v2).y * (v3).y, (v1).z = (v2).z * (v3).z
#define VMulTwo(v1, v2, v3)  (v1).x = (v2).x * (v3).x, (v1).y = (v2).y * (v3).y
#define VDiv(v1, v2, v3)  (v1).x = (v2).x / (v3).x, (v1).y = (v2).y / (v3).y, (v1).z = (v2).z / (v3).z
#define VDivTwo(v1, v2, v3)  (v1).x = (v2).x / (v3).x, (v1).y = (v2).y / (v3).y  

// Sum (square v1 and divide by v2, element by element) - useful to compute rotational kinEnergy: 
#define VSumSqDiv(v1, v2)  ((v1).x * (v1).x / (v2).x + (v1).y * (v1).y / (v2).y + (v1).z * (v1).z / (v2).z) 

#define VSAdd(v1, v2, s3, v3)  (v1).x = (v2).x + (s3) * (v3).x, (v1).y = (v2).y + (s3) * (v3).y, (v1).z = (v2).z + (s3) * (v3).z
#define VSAddTwo(v1, v2, s3, v3)  (v1).x = (v2).x + (s3) * (v3).x, (v1).y = (v2).y + (s3) * (v3).y
#define VSSAdd(v1, s2, v2, s3, v3) (v1).x=(s2)*(v2).x+(s3)*(v3).x,(v1).y=(s2)*(v2).y+(s3)*(v3).y,(v1).z=(s2)*(v2).z+(s3)*(v3).z

/* product operations: */
#define VDot(v1, v2)  ((v1).x * (v2).x + (v1).y * (v2).y + (v1).z * (v2).z)
#define VDotTwo( v1, v2 )  ( ( v1 ).x * ( v2 ).x + ( v1 ).y * ( v2 ).y )
#define VWDot(v1, v2, v3)  ((v1).x * (v2).x * (v3).x + (v1).y * (v2).y * (v3).y + (v1).z * (v2).z * (v3).z)
#define VCross(v1, v2, v3)  (v1).x = (v2).y * (v3).z - (v2).z * (v3).y,	(v1).y = (v2).z * (v3).x - (v2).x * (v3).z, (v1).z = (v2).x * (v3).y - (v2).y * (v3).x

/* matrix-vector products (for both the matrix and its transpose): */
#define MVMul(v1, m, v2)						\
  (v1).x = (m)[0] * (v2).x + (m)[3] * (v2).y + (m)[6] * (v2).z,		\
    (v1).y = (m)[1] * (v2).x + (m)[4] * (v2).y + (m)[7] * (v2).z,	\
    (v1).z = (m)[2] * (v2).x + (m)[5] * (v2).y + (m)[8] * (v2).z
#define MVMulT(v1, m, v2)						\
  (v1).x = (m)[0] * (v2).x + (m)[1] * (v2).y + (m)[2] * (v2).z,		\
    (v1).y = (m)[3] * (v2).x + (m)[4] * (v2).y + (m)[5] * (v2).z,	\
    (v1).z = (m)[6] * (v2).x + (m)[7] * (v2).y + (m)[8] * (v2).z

/* other useful operations */
#define MCopy( m1, m2 )						\
  (m1)[0] = (m2)[0], (m1)[1] = (m2)[1], (m1)[2] = (m2)[2],	\
    (m1)[3] = (m2)[3], (m1)[4] = (m2)[4], (m1)[5] = (m2)[5],	\
    (m1)[6] = (m2)[6], (m1)[7] = (m2)[7], (m1)[8] = (m2)[8]

#define MScale(m, s)				\
  (m)[0] *= s,(m)[3] *= s,(m)[6] *= s,		\
    (m)[1] *= s,(m)[4] *= s,(m)[7] *= s,	\
    (m)[2] *= s,(m)[5] *= s,(m)[8] *= s

/* other useful operations */
#define MExpScale(m1, s, m2)						\
  (m1)[0] *= exp(s * (m2)[0]), (m1)[3] *= exp(s * (m2)[3]), (m1)[6] *= exp(s * (m2)[6]), \
    (m1)[1] *= exp(s * (m2)[1]), (m1)[4] *= exp(s * (m2)[4]), (m1)[7] *= exp(s * (m2)[7]), \
    (m1)[2] *= exp(s * (m2)[2]), (m1)[5] *= exp(s * (m2)[5]), (m1)[8] *= exp(s * (m2)[8])

// Matrix Determinant, assuming a 3 x 3 matrix
#define MatDet( m )							\
  (m)[0]*(m)[4]*(m)[8] - (m)[0]*(m)[7]*(m)[5] - (m)[3]*(m)[1]*(m)[8] + (m)[3]*(m)[7]*(m)[2] + (m)[6]*(m)[1]*(m)[5] - (m)[6]*(m)[4]*(m)[2]  

// Matrix addition - element by element
#define MAdd( m1, m2, m3 )			\
  (m1)[0] = (m2)[0] + (m3)[0],			\
    (m1)[1] = (m2)[1] + (m3)[1],		\
    (m1)[2] = (m2)[2] + (m3)[2],		\
    (m1)[3] = (m2)[3] + (m3)[3],		\
    (m1)[4] = (m2)[4] + (m3)[4],		\
    (m1)[5] = (m2)[5] + (m3)[5],		\
    (m1)[6] = (m2)[6] + (m3)[6],		\
    (m1)[7] = (m2)[7] + (m3)[7],		\
    (m1)[8] = (m2)[8] + (m3)[8]

// Matrix subtraction - element by element
#define MSub( m1, m2, m3 )			\
  (m1)[0] = (m2)[0] - (m3)[0],			\
    (m1)[1] = (m2)[1] - (m3)[1],		\
    (m1)[2] = (m2)[2] - (m3)[2],		\
    (m1)[3] = (m2)[3] - (m3)[3],		\
    (m1)[4] = (m2)[4] - (m3)[4],		\
    (m1)[5] = (m2)[5] - (m3)[5],		\
    (m1)[6] = (m2)[6] - (m3)[6],		\
    (m1)[7] = (m2)[7] - (m3)[7],		\
    (m1)[8] = (m2)[8] - (m3)[8]

// Matrix accumulation sum: m1 += m2 (element by element)
#define MMAdd( m1, m2 )				\
  (m1)[0] += (m2)[0],				\
    (m1)[1] += (m2)[1],				\
    (m1)[2] += (m2)[2],				\
    (m1)[3] += (m2)[3],				\
    (m1)[4] += (m2)[4],				\
    (m1)[5] += (m2)[5],				\
    (m1)[6] += (m2)[6],				\
    (m1)[7] += (m2)[7],				\
    (m1)[8] += (m2)[8]

// Matrix accumulation subtraction: m1 -= m2 (element by element)
#define MMSub( m1, m2 )				\
  (m1)[0] -= (m2)[0],				\
    (m1)[1] -= (m2)[1],				\
    (m1)[2] -= (m2)[2],				\
    (m1)[3] -= (m2)[3],				\
    (m1)[4] -= (m2)[4],				\
    (m1)[5] -= (m2)[5],				\
    (m1)[6] -= (m2)[6],				\
    (m1)[7] -= (m2)[7],				\
    (m1)[8] -= (m2)[8]

// Matrix addition with scaling - element by element
#define MSAdd( m1, m2, s3, m3 )						\
  (m1)[0] = (m2)[0] + (s3)*(m3)[0],  (m1)[3] = (m2)[3] + (s3)*(m3)[3],  (m1)[6] = (m2)[6] + (s3)*(m3)[6], \
    (m1)[1] = (m2)[1] + (s3)*(m3)[1],  (m1)[4] = (m2)[4] + (s3)*(m3)[4],  (m1)[7] = (m2)[7] + (s3)*(m3)[7], \
    (m1)[2] = (m2)[2] + (s3)*(m3)[2],  (m1)[5] = (m2)[5] + (s3)*(m3)[5],  (m1)[8] = (m2)[8] + (s3)*(m3)[8]

#define MMSAdd( m1, s2, m2 )  MSAdd( m1, m1, s2, m2 )

#define MZero( m ) (m)[0] = 0, (m)[1] = 0, (m)[2] = 0, (m)[3] = 0, (m)[4] = 0, (m)[5] = 0, (m)[6] = 0, (m)[7] = 0, (m)[8] = 0

// Set off-diagonal elements to zero
#define MZeroOffDiags( m )   (m)[1] = 0, (m)[2] = 0, (m)[3] = 0, (m)[5] = 0, (m)[6] = 0, (m)[7] = 0

// Trace operator
#define Trace( m ) ( (m)[0] +  (m)[4] +  (m)[8] )

// Outer product between vectors v1, v2 - result is (stored in) matrix m 
#define VOuter( m, v1, v2 )			\
  (m)[0] = (v1).x * (v2).x,			\
    (m)[1] = (v1).y * (v2).x,			\
    (m)[2] = (v1).z * (v2).x,			\
    (m)[3] = (v1).x * (v2).y,			\
    (m)[4] = (v1).y * (v2).y,			\
    (m)[5] = (v1).z * (v2).y,			\
    (m)[6] = (v1).x * (v2).z,			\
    (m)[7] = (v1).y * (v2).z,			\
    (m)[8] = (v1).z * (v2).z  

#define VProd(v)				\
  ((v).x * (v).y * (v).z)
#define VGe(v1, v2)						\
  ((v1).x >= (v2).x && (v1).y >= (v2).y && (v1).z >= (v2).z)
#define VLt(v1, v2)						\
  ((v1).x < (v2).x && (v1).y < (v2).y && (v1).z < (v2).z)
#define VLinear(p, s)				\
  (((p).z * (s).y + (p).y) * (s).x + (p).x)
#define VSetAll(v, s)				\
  VSet (v, s, s, s)
#define VSetAllTwo(v, s)   VSetTwo (v, s, s)
#define VAddCon(v1, v2, s)			\
  (v1).x = (v2).x + (s),			\
    (v1).y = (v2).y + (s),			\
    (v1).z = (v2).z + (s)
#define VCSum(v)				\
  ((v).x + (v).y + (v).z)
#define VComp(v, k)					\
  *((k == 0) ? &(v).x : ((k == 1) ? &(v).y : &(v).z))

/* two operations for converting between vectors and array elements */
#define VToLin(a, n, v)				\
  a[(n) + 0] = (v).x,				\
    a[(n) + 1] = (v).y,				\
    a[(n) + 2] = (v).z
#define VFromLin(v, a, n)			\
  VSet (v, a[(n) + 0], a[(n) + 1], a[(n) + 2])


/* additional vector operations */
#define VZeroTwo(v)  VSetAllTwo (v, 0)
#define VZero(v)  VSetAll (v, 0)
#define VLenSq(v)  VDot (v, v)
#define VLenSqTwo(v)  VDotTwo( v, v )
#define VWLenSq(v1, v2)  VWDot(v1, v2, v2)
#define VLen(v)  sqrt (VDot (v, v))
#define VVAdd(v1, v2)  VAdd (v1, v1, v2)
#define VVAddTwo(v1, v2)  VAddTwo (v1, v1, v2)
#define VVSub(v1, v2)  VSub (v1, v1, v2)
#define VVSAdd(v1, s2, v2) VSAdd (v1, v1, s2, v2)
#define VVSAddTwo(v1, s2, v2) VSAddTwo (v1, v1, s2, v2)
#define VInterp(v1, s2, v2, v3)			\
  VSSAdd (v1, s2, v2, 1. - (s2), v3)

/* quaternions: definition of dataStructs and operations */
typedef struct {
  real u1, u2, u3, u4;
} Quat;

#define QSet(q, s1, s2, s3, s4)			\
  (q).u1 = s1,					\
    (q).u2 = s2,				\
    (q).u3 = s3,				\
    (q).u4 = s4
#define QZero(q)  QSet (q, 0, 0, 0, 0)
#define QScale(q, s)				\
  (q).u1 *= s,					\
    (q).u2 *= s,				\
    (q).u3 *= s,				\
    (q).u4 *= s
#define QSAdd(q1, q2, s3, q3)			\
  (q1).u1 = (q2).u1 + (s3) * (q3).u1,		\
    (q1).u2 = (q2).u2 + (s3) * (q3).u2,		\
    (q1).u3 = (q2).u3 + (s3) * (q3).u3,		\
    (q1).u4 = (q2).u4 + (s3) * (q3).u4
#define QLenSq(q)				\
  (Sqr ((q).u1) + Sqr ((q).u2) + Sqr ((q).u3) +	\
   Sqr ((q).u4))
#define QMul(q1, q2, q3)				\
  (q1).u1 =   (q2).u4 * (q3).u1 - (q2).u3 * (q3).u2 +	\
    (q2).u2 * (q3).u3 + (q2).u1 * (q3).u4,		\
    (q1).u2 =   (q2).u3 * (q3).u1 + (q2).u4 * (q3).u2 -	\
    (q2).u1 * (q3).u3 + (q2).u2 * (q3).u4,		\
    (q1).u3 = - (q2).u2 * (q3).u1 + (q2).u1 * (q3).u2 +	\
    (q2).u4 * (q3).u3 + (q2).u3 * (q3).u4,		\
    (q1).u4 = - (q2).u1 * (q3).u1 - (q2).u2 * (q3).u2 -	\
    (q2).u3 * (q3).u3 + (q2).u4 * (q3).u4

/* complex variables */
typedef struct {
  real R, I;
} Cmplx;

#define CSet(a, x, y)				\
  a.R = x,					\
    a.I = y
#define CAdd(a, b, c)				\
  a.R = b.R + c.R,				\
    a.I = b.I + c.I
#define CSub(a, b, c)				\
  a.R = b.R - c.R,				\
    a.I = b.I - c.I
#define CMul(a, b, c)				\
  a.R = b.R * c.R - b.I * c.I,			\
    a.I = b.R * c.I + b.I * c.R

/* From Rapport, [p.498] - constant used in PutConfig() note: 2^(15)=32768... */
#define SCALE_FAC 32767 

/* macro introduced to aid readability by hiding the fact that matrices
   are in fact stored as singly-indexed array */
#define MAT( a, n, i, j )  ( a )[ ( i ) + n * ( j ) ]

/* macros handling allocations of one- and two-dimensional arrays of
   any kind of variable or structure */
#define AllocMem(a, n, t)  a = (t *) malloc ((n) * sizeof (t))
#define AllocMem2(a, n1, n2, t)				\
  AllocMem (a, n1, t *);				\
  AllocMem (a[0], (n1) * (n2), t);			\
  for (k = 1; k < n1; k ++) a[k] = a[k - 1] + n2;

#define MAX_MPEX_ORD  2
#define I(i, j)  ((i) * ((i) + 1) / 2 + (j))
#define c(i, j)  c[I(i, j)]
#define s(i, j)  s[I(i, j)]

typedef struct {
  real c[I(MAX_MPEX_ORD, MAX_MPEX_ORD) + 1], s[I(MAX_MPEX_ORD, MAX_MPEX_ORD) + 1];
} MpTerms;
typedef struct {
  MpTerms le, me;
  int occ;
} MpCell;

/* shorthand for frequently appearing loops */
#define DO_SITE          for ( n = 0; n < nSites; n++ )
#define DO_ATOM          for ( a = 0; a < nAtoms; a++ )
#define DO_WATER         for ( n = 0; n < nWaters; n++ )
#define DO_LIPID         for ( n = 0; n < nLipids; n++ )
#define DO_DMPCLIPID      for ( s = 0; s < dmpcStruct.nBeads; s++ ) 
#define DO_CELL( j, m )  for ( j = cellList[ m ]; j >= 0; j = cellList[ j ] )

/* macros dealing with Periodic Boundary Conditions - TO BE MODIFIED FOR NON-ISOSHAPE NPT DYNAMICS!!! */
#define VWrap(v, t)					\
  if (v.t >= 0.5 * region.t)      v.t -= region.t;	\
  else if (v.t < -0.5 * region.t) v.t += region.t

#define VWrapS(v, t)				\
  if (v.t >= 0.5 )      v.t -= 1.;		\
  else if (v.t < -0.5 ) v.t += 1.

#define VShift(v, t)					\
  if (v.t >= 0.5 * region.t)      shift.t -= region.t;	\
  else if (v.t < -0.5 * region.t) shift.t += region.t

#define VShiftWrap(v, t)			\
  if (v.t >= 0.5 * region.t) {			\
    shift.t -= region.t;			\
    v.t -= region.t;				\
  } else if (v.t < -0.5 * region.t) {		\
    shift.t += region.t;			\
    v.t += region.t;				\
  }

#define VCellWrap(t)				\
  if (m2v.t >= cells.t) {			\
    m2v.t = 0;					\
    shift.t = region.t;				\
  } else if (m2v.t < 0) {			\
    m2v.t = cells.t - 1;			\
    shift.t = - region.t;			\
  }

#define VWrapAll(v)				\
  {VWrap (v, x);				\
    VWrap (v, y);				\
    VWrap (v, z);}

#define VWrapXY( v ) { VWrap( v, x );  VWrap( v, y ); }

#define VWrapAllS(v)				\
  {VWrapS (v, x);				\
    VWrapS (v, y);				\
    VWrapS (v, z);}
#define VShiftAll(v)				\
  {VShift (v, x);				\
    VShift (v, y);				\
    VShift (v, z);}
#define VCellWrapAll()				\
  {VCellWrap (x);				\
    VCellWrap (y);				\
    VCellWrap (z);}

/* cell offsets used for interaction and neighbour-list computations */
#define OFFSET_VALS							\
  { {0,0,0}, {1,0,0}, {1,1,0}, {0,1,0}, {-1,1,0},			\
					  {0,0,1}, {1,0,1}, {1,1,1}, {0,1,1}, {-1,1,1}, {-1,0,1}, \
											  {-1,-1,1}, {0,-1,1}, {1,-1,1}	\
  }

#define N_OFFSET  14

/* measurements stuff */
typedef struct { /* property measurements */
  real val, /* actual measured value */ 
    sum,    /* sum accumulated, when averaged gives the mean */
    sum2;   /* sum of squares, when averaged gives the standard deviation */
} Prop;

/*********************************************************************************************
 * The following operations are defined to be used with the "Prop" structures [Rapaport, p.32]
 *********************************************************************************************/
// PropZero initialises
#define PropZero(v)  v.sum = v.sum2 = 0. 
// PropAccum accumulates
#define PropAccum(v)  v.sum += v.val, v.sum2 += Sqr (v.val)
// PropAvg calculates the mean (v.sum /= n) and v.sum2 = sqrt( variance ) = standard deviation = rms fluctuation
#define PropAvg(v, n)  v.sum /= n, v.sum2 = sqrt (Max (v.sum2 / n - Sqr (v.sum), 0.))
// PropEst: property estimate in the format (mean, standard deviation) 
#define PropEst(v)  v.sum, v.sum2
// PropScaleEst: similar to PropEst, but applies scaling factor - useful to convert output to physical units  
#define PropScaleEst( v, s )  v.sum * s, v.sum2 * s 
// PropShiftEst: similar to PropEst, but applies shifting factor - useful to convert output to physical units  
#define PropShiftEst( v, s )  v.sum + s, v.sum2  

#define SITE_LIMIT  10000000

#define NameVal(x)				\
  if (! strncmp (bp, #x, strlen (#x))) {	\
    bp += strlen (#x);				\
    x = strtod (bp, &bp);			\
  }

#define NHIST  (5)  /* NDIM + 2 */

/* file stuff */
#define ReadF( x )       fread( &x, sizeof( x ), 1, fp )
#define WriteF( x )      fwrite( &x, sizeof( x ), 1, fp )
#define ReadFN( x, n )   fread( x, sizeof( x[ 0 ] ), n, fp )
#define WriteFN( x, n )  fwrite( x, sizeof( x[ 0 ] ), n, fp )

enum errCodes {ERR_NONE, ERR_BOND_SNAPPED, ERR_CHECKPT_READ, ERR_CHECKPT_WRITE,
	       ERR_COPY_BUFF_FULL, ERR_EMPTY_EVPOOL, ERR_MSG_BUFF_FULL,
	       ERR_OUTSIDE_REGION, ERR_SNAP_READ, ERR_SNAP_WRITE,
	       ERR_SUBDIV_UNFIN, ERR_TOO_MANY_CELLS, ERR_TOO_MANY_COPIES,
	       ERR_TOO_MANY_LAYERS, ERR_TOO_MANY_LEVELS, ERR_TOO_MANY_SITES,
	       ERR_TOO_MANY_MOVES, ERR_TOO_MANY_NEBRS, ERR_TOO_MANY_REPLICAS, 
               ERR_UNDETERMINED_INTERACTION, ERR_STRUCT_BUILT_INCONSISTENT_WITH_INPUT};

#endif
